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Abstract. The formation of local denaturation zones (bubbles) in double-stranded 
DNA is an important example for conformational changes of biological macromolecules. 
We study the dynamics of bubble formation in terms of a Fokker-Planck equation for 
the probability density to find a bubble of size n base pairs at time t, on the basis of the 
free energy in the Poland-Scheraga model. Characteristic bubble closing and opening 
times can be determined from the corresponding first passage time problem, and are 
sensitive to the specific parameters entering the model. A multistate unzipping model 
with constant rates recently applied to DNA breathing dynamics [G. Altan-Bonnet et 
al , Phys. Rev. Lett. 90, 138101 (2003)] emerges as a limiting case. 
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Introduction. Under physiological conditions, the equilibrium structure of a DNA 
molecule is the double-stranded Watson-Crick helix. At the same time, in essentially 
all physiological processes involving DNA, for docking to the DNA, DNA binding 
proteins require access to the "inside" of the double-helix, and therefore the unzipping 
(denaturation) of a specific region of base pairs [HE]- Examples include the replication 
of DNA via DNA helicase and polymerase, and transcription to single-stranded DNA 
via RNA polymerase. Thus, double-stranded DNA has to open up locally to expose the 
otherwise satisfied bonds between complementary bases. 

There are several mechanisms how such unzipping of double-stranded DNA 
(dsDNA) can be accomplished. Under physiological conditions, local unzipping occurs 
spontaneously due to fluctuations, the breathing of dsDNA, which opens up bubbles 
of a few tens of base pairs These breathing fluctuations may be supported by 
accessory proteins which bind to transient single-stranded regions, thereby lowering the 
DNA base pair stability pj. Single molecule force spectroscopy opens the possibility 
to induce denaturation regions of controllable size, by pulling the DNA with optical 
tweezers [IJ. In this way, the destabilising activity of the ssDNA-binding T4 gene 32 
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protein has been probed, and a kinetic barrier for the single-strand binders identified 
Denaturation bubbles can also be induced by under-winding the DNA double helix 
[Bj. A recent study of the dynamics of these twist-induced bubbles in a random DNA 
sequence shows that small bubbles (less than several tens of base pairs) are delocalised 
along the DNA, whereas larger bubbles become localised in AT-rich regions [Zj. Finally, 
upon heating, dsDNA exhibits denaturation bubbles of increasing size and number, and 
eventually the two strands separate altogether in a process called denaturation transition 
or melting [HI E]. Depending on the specific sequence of the DNA molecule and the 
solvent conditions, the temperature T m at which one-half of the DNA has denatured 
typically ranges between 50° C and 100° C. Due to the thermal lability of typical natural 
proteins, thermal melting of DNA is less suited for the study of protein-DNA interactions 
than force-induced melting 4j. On the other hand, the controlled melting of DNA by 
heating is an important step of the PCR method for amplifying DNA samples ^0], with 
numerous applications in biotechnology [TT] . 

The study of the bubble dynamics in the above processes is of interest in view of 
a better understanding of the interaction with single-stranded DNA binding proteins. 
This interaction involves an interplay between different time scales, e.g., the relaxation 
time of the bubbles and the time needed for the proteins to rearrange sterically in 
order to bind ^2]- Dynamic probes such as single-molecule force spectroscopy jlj and 
molecular beacon assays ^Hj may therefore shed light on the underlying biochemistry 
of such processes. 

In a recent experiment by Altan-Bonnet et al ^1] , the dynamics of a single bubble 
in dsDNA was measured by fluorescence correlation spectroscopy. It was found that in 
the breathing domain of the DNA construct (a row of 18 AT base pairs sealed by more 
stable GC base pairs) fluctuation bubbles of size 2 to 10 base pairs are formed below the 
melting temperature T m of the AT breathing domain. The relaxation dynamics follows 
a multistate relaxation kinetics involving a distribution of bubble sizes and successive 
opening and closing of base pairs. The characteristic relaxation time scales were 
estimated from the experiment to within the range of 20 to 100 fis. Also in reference |14j . 
a simple master equation of stepwise zipping-unzipping with constant rate coefficients 
was proposed to successfully describe the data for the autocorrelation function, showing 
that indeed the bubble dynamics is a multistate process. The latter was confirmed in a 
recent UV light absorption study of the denaturation of DNA oligomers ^3] . 

In this work, we establish a general framework to study the bubble dynamics of 
dsDNA by means of a Fokker-Planck equation, based on the bubble free energy function. 
The latter allows one to include microscopic interactions in a straightforward fashion, 
such that our approach may serve as a testing ground for different models, as we show 
below. In particular, it turns out that the phenomenological rate equation approach, 
corresponding to a diffusion with constant drift in the space of bubble size n used 
by Altan-Bonnet et al to fit their experimental data corresponds to a limiting 
case of our Fokker-Planck equation. However, the inclusion of additional microscopic 
interactions in such a rate equation approach is not straightforward I n what follows, 
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we first establish the bubble free energy within the Poland- Scheraga (PS) model of DNA 
denaturation [HI El; and then derive the Fokker-Planck equation to describe the bubble 
dynamics both below and at the melting temperature of dsDNA. 

Bubble free energy. In the PS model, energetic bonds in the double-stranded, helical 
regions of the DNA compete with the entropy gain from the far more flexible, single- 
stranded loops [SI E] • The stability of the double helix originates mainly from stacking 
interactions between adjacent base pairs, aside from the Watson-Crick hydrogen bonds 
between bases. In addition, the positioning of bases for pairing out of a loop state 
gives rise to an entropic contribution. The Gibbs free energy Gjj = Hij — TSij for the 
dissociation of two paired and stacked base pairs i and j has been measured, and is 
available in terms of the enthalpic and entropic contributions and SV, 17 . In the 
following we consider a homopolymer for simplicity, as suitable for the AT breathing 
domain in reference [Hj. For an AT- homopolymer (i = j = AT), the Gibbs free 
energy per base pair in units of fcgT yields as 7 = j3Gu/2 = 0.6 at 37°C for standard 
salt conditions (0.0745 M-Na + ). Similarly, for a GC- homopolymer one finds the higher 
value of 7 = 1.46 at 37°C. The condition 7 = defines the melting temperature T m 
HZIIig, thus T m (AT) = 66.8°C and T m (GC) = 102.5°C (we assume that is linear in 
T, cf. reference jT9*]). Above T m , 7 becomes negative. For given 7 = 7(T), the statistical 
weight for the dissociation of n base pairs obtains as 



Additional contributions arise upon formation of a loop within dsDNA. Firstly, an initial 
energy barrier has to be overcome, which we denote as 71 in units of ksT. From fitting 
melting curves to long DNA, 71 « 10 was obtained, so that the statistical weight for the 
initiation of a loop (cooperativity parameter), a = exp(— 71), is of order 10~ 5 jSlEj- As 
the energy to extend an existing loop by one base pair is smaller than fceT, DNA melts 
as large cooperative domains. Below the melting temperature T m , the bubbles become 
smaller, and long range interactions beyond nearest neighbours become more important. 
In this case, the probability of bubble formation is larger, and 71 ranges between 3 and 
5, thus cr ;$ 0.05 [7] (cf. section 5 in [9 ). According to reference [TH] , the smallness of 
(Jo inhibits the recombination of complementary DNA strands with mutations, making 
recognition more selective. Secondly, once a loop of n base pairs has formed, there is a 
weight f(n) of mainly entropic origin, to be detailed below. The additional weight of a 
loop of n open base pairs is thus 



For large bubbles one usually assumes the form /(n) = (n+ l)~ c 0HII. Here, the value 
of the exponent c = 1.76 corresponds to a self-avoiding, flexible ring [HI EH ED] , which 
is classically used in denaturation modelling within the PS approach. Recently, the 
PS model has been considered in view of the order of the denaturation phase transition 
[2U E2J EHJ EU EE| • Reference ^H] finds by finite size scaling analysis of measured melting 
curves of DNA oligomers that the transition is of second order. In reference [2J, the 
value c = 2.12 was suggested, compare the discussion in references 0EUEH]- For smaller 



W(n) = exp(— wy). 



(1) 



fi(w) = o- f(n) . 



(2) 
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bubbles (in the range of 1 to a few tens of base pairs), the appropriate form of fin) is 
more involved. In particular, f(n) will depend on the finite persistence length of ssDNA 
(about eight bases), on the specific base sequence, and possibly on interactions between 
dissolved but close-by base pairs (cf. section 2.1.3 in 9J). Therefore, the knowledge of 
f(n) provides information on these microscopic interactions. 

For simplicity, we here adopt the simple form f(n) — (n + l)~ c for all n > 0, and 
consider the loop weight 

Q(n) = a (n + l)- c . (3) 

We show that at the melting temperature the results for the relaxation times for the 
bubbles are different for the available values c = 1.76 and c = 2.12 quoted above. This 
shows that the specific form of f(n) indeed modifies experimentally accessible features 
of the bubble dynamics 14 . 

In what follows, we focus on a single bubble in dsDNA, neglecting its interaction 
with other bubbles. Since due to er <C 1 the mean distance between bubbles (~ l/cr 
[Hj) is large, this approximation is justified as long as the bubbles are not too large [Zj. 
It also corresponds to the situation studied in the recent experiment by Altan-Bonnet 
et al [H]. According to the above, the total free energy Fin) of a single bubble with 
n > open base pairs follows in the form 

PF(n) = - In [W(n)Q(n)] = nj(T) + 7l + c\n(n + 1) , (4) 

where the dependence on the temperature T enters only via 7 = 7 (T). We show the 
free energy (p|J) in figure Q for c = 1.76 and for the parameters of an AT-homopolymer, 
for physiological temperature T = 37°C, at the melting temperature T m = 66°C, and 
at T = 100°, compare the discussion below. 

Bubble dynamics. In the generally accepted multistate unzipping model, the double 
strand opens by successively breaking Watson-Crick bonds, like opening a zipper [2l)|l27j. 
As 7 becomes small on increasing the temperature, thermal fluctuations become relevant 
and cause a random walk-like propagation of the zipper locations at both ends of the 
bubble-helix joints. The fluctuations of the bubble size can therefore be described 
in the continuum limit through a Fokker-Planck equation for the probability density 
function (PDF) P(n,t) to find at time t a bubble consisting of n denatured base pairs, 
following a similar reasoning as pursued in the modelling of the dynamics of biopolymer 
translocation through a narrow membrane pore To establish this Fokker-Planck 

equation, we combine the continuity equation (compare reference |28| ) 

dP(n,t) dj(n,t) 

-^z- + ^n-=° (5) 
with the expression for the corresponding flux, 

where it is assumed that the potential exerting the drift is given by the bubble free energy 
(J1J. In equation (jSJ, we incorporated an Einstein relation of the form D = ksTfj,, where 
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Figure 1. Variable part of the bubble free energy Q, (3J-{n) — 71 = n-f(T)+c\n(n+l), 
as a function of the bubble size n for T = 37°C (7 w 0.6), T m = 66°C (7 = 0), and 
T = 100°C (7 w -0.5). In the latter small barrier precedes the negative drift 

towards bubble opening as dominated by the 7 < contribution. 



the mobility fi has dimensions \p] = sec/(g • cm 2 ), and therefore [D] = sec -1 represents 
an inverse time scale. By combination of equations (J1J, © and (|5|l. we retrieve the 
Fokker-Planck equation for P(n,t) f 



Thus, we arrived at a reduced ID description of the bubble dynamics in a homopolymer, 
with the bubble size n as the effective 'reaction' coordinate. For a heteropolymer, also 
the position of the bubble within the double helix, i.e., the index m of the first open base 
pair, becomes relevant. In this case, the bubble free energy T and thus the PDF depend 
on both m and n, and on the specific base sequence; the corresponding generalisation 
of equation (JJJ) is straightforward. In a random sequence, additional phenomena may 
occur, such as localisation of larger bubbles Finally, to establish the Fokker-Planck 
equation (J7j), we assume that changes of the bubble size n occur slower than other 
degrees of freedom of the PS free energy within the bubble region (e.g., Rouse-Zimm 
modes). This assumption seems reasonable due to the long bubble dynamics' relaxation 
time scales of 20 to 100 /is [14J, and the good approximation of bubble independence [7j. 

By rescaling time according to t — » Dt, the Fokker-Planck equation (JJJ) can be 
made dimensionless, a representation we are going to use in the numerical evaluation 
below. The formulation in terms of a Fokker-Planck equation makes it possible to derive 
the characteristic times for bubble closing and opening in terms of a first passage time 
problem. That is, for bubble closing, the associated mean closing time follows as the 
mean first passage time to reach bubble size n = after starting from the initial bubble 
size no- We now determine these characteristic times for the three regimes defined by 7 
with respect to temperature T. 

I Note that the operator acts also on P(n,t). 
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Figure 2. Characteristic bubble closing times r as a function of initial bubble size no 
for an AT-homopolymer, obtained from the Fokker-Planck equation Q by numerical 
integration. At T — 37°C, the result for c = 1.76 is compared to the approximation 
c = which leads to somewhat larger closing times. The analytical solution for 
c = compares well with the numerical result, the slight discrepancy being due to the 
reflecting boundary condition applied in the numerics, in comparison to the natural 
boundary condition at n — > oo used to derive equation (JHJ) . At the melting temperature 
T m = 66°C, the closing times for the values c = 1.76 and c = 2.12 can be distinguished. 



(i) T < T m . In this regime, the drift consists of two contributions, the constant 
drift -D7 and the loop closure component Dc/(n + 1), which decreases with n. For 
large n, we can therefore approximate the drift by the constant term Z>y, and in this 
limit the Fokker-Planck equation (jZJ) is equivalent to the continuum limit of the master 
equation used in reference ^3] to describe the experimental bubble data. In particular, 
we can identify our two independent parameters D and 7 with the rate constants k + 
and to open and close a base pair introduced in jHj, respectively: D = {k + + fc_)/2 
and 7 = 2{k_ — k + )/(k + + kJ). In this approximation, the correlation functions used 
successfully to fit the experimental results in ^1] can be derived from the Fokker-Planck 
equation (|7jl. Moreover, we can deduce the mean first passage time PDF for a bubble 
of initial size no to close in the exact analytical form (compare |29|) 

V ; V^LW 1 ^Dt J v ; 

which decays exponentially for large n. In particular, from (JBJ) the characteristic (mean) 
first passage time for bubble closing, 

r = n /(£> 7 ) (9) 

follows, which is linear in the initial bubble size n . In figure 121 we compare this 
analytical result for the value 7(37°C) with the characteristic closing times using the 
full drift term from equation (JJJ), obtained from numerical integration. It can be seen 
that the qualitative behaviour for both cases with and without the Dc/ (n + 1) term is 
very similar, but that in the presence of the loop closure component the characteristic 
times are reduced. 
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(ii.) T = T m . At the melting temperature, the drift in equation (J7J) is solely given 
by the loop closure term Dc/(n + 1). In figure 121 we plot the characteristic bubble 
closing time obtained numerically. In comparison to the above case T < T m , the faster 
than linear increase as function of initial bubble size n is distinct. Keeping in mind 
that Dcj (n+ 1) becomes very small for increasing n, this behaviour can be qualitatively 
understood from the approximation in terms of a drift-free diffusion in a box of size n , 
in which the initial condition P(n, 0) = 5(n — n ) is located at a reflecting boundary, 
and at n = an absorbing boundary is placed. This problem can be solved analytically, 
with the result r = n1/{2D) for the characteristic escape time § in which the quadratic 
dependence on uq contrasts the linear behaviour in the result 0. Thus, at the melting 
temperature T m , the tendency for a bubble to close becomes increasingly weaker for 
larger bubble sizes, and therefore much larger bubbles can be formed, in contrast to the 
case T < T m . A further comparison to the value c = 2.12 mentioned above demonstrates 
that a clear quantitative difference in the associated closing times exists, see figure El 
However, the qualitative behaviour remains unchanged. In principle, the study of the 
bubble dynamics can therefore be used to discern different models for the loop closure 
factor. 

(iii.) T > T m . Above the melting temperature, the drift is governed by the interplay 
between the loop closure component Dcj (n + 1) tending to close the bubble, and the 
bubble free energy D^{T > T m ) < 0, which causes a bias towards bubble opening. In 
figure ^ we show for the AT-homopolymer case how the overall drift potential after 
a small initial activation barrier becomes negative, and the dynamics is essentially 
governed by the 7-contribution. As a consequence, from the result (jHJ), we find that 
the associated mean first passage time diverges, i.e., the bubble on average increases in 
size until the entire DNA is denatured, as expected from the PS model. 

By symmetry, similar results hold for the bubble opening process. However, the 
existence of the bubble initiation energy 71 involves an additional Arrhenius factor, 
which is not included in the Fokker-Planck equation (J7|), and which reduces the opening 
probability, causing an increase of the associated opening time, cf. reference |14j . 

Conclusions. From the bubble free energy of a single, independent bubble in the 
Poland- Scheraga theory of DNA melting, we derived a Fokker-Planck equation for the 
PDF to find a bubble created by denaturation of n base pairs at time t. This formulation 
allows for the calculation of the characteristic times scales for bubble closing and opening 
in terms of first passage time problems. Three different regimes were distinguished: (i) 
below the melting temperature, the characteristic bubble closing time increases linearly 
in bubble size, and the drift can be approximated by the constant value -D7. In this 
approximation, the Fokker-Planck equation matches the continuum version of the master 
equation employed previously in an experimental study of DNA breathing |14j . (ii) At 
the melting temperature, an approximately quadratic growth of the bubble closing time 
as a function of bubble size is observed, which can be explained by noting that the 

§ This result can, for instance, be obtained through the method of images, compare reference 
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l/(n + Independence of the loop closure drift component can be neglected for larger 
n, leading to pure diffusion. In this approximation, the exact analytical results indeed 
lead to the quadratic dependence observed in the numerical results, (iii) Above the 
melting temperature, the characteristic closing time in our model diverges, consistent 
with the fact that on average the DNA follows the trend towards the thermodynamically 
favourable state of complete denaturation. 

The expression for the Gibbs free energy used in our approach involves a purely 
entropic contribution for a single-stranded bubble. It was suggested in reference ^1] 
that also in denaturation bubbles a residual stacking energy e s is present. In this case, 
the value of 7 would have to be corrected by this e s . 

From the Fokker-Planck equation (|7J), in which the microscopic interactions enter 
via the free energy one can derive measurable quantities such as moments and 
dynamic correlation functions In principle, the Fokker-Planck equation involves one 
free parameter, the time scale 1/D, while the values for the other parameters are known. 
However, by fitting sufficiently accurate experimental data for DNA bubble dynamics at 
different temperatures, the values for additional parameters may be extracted by help 
of our general Fokker-Planck framework probing suitable free energy functions. 
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